Subsets only the NK and T cells.
library(Seurat)
## Attaching SeuratObject
## Attaching sp
library(dittoSeq)
## Loading required package: ggplot2
library(readxl)
sc <- readRDS("analysis/sc_integrated_annotated.rds")
tc <- subset(sc, subset = celltype_fine %in% c("NK cells", "T cells"))
rm(sc)
# load in file that can be used to re-orde samples in e.g. barplot
sample_order <- read_excel("metadata/sample_order.xlsx")
sample_order$sample <- sapply(sample_order$`Sample ID`, function(x){
strsplit(x, " ")[[1]][1]
})
DefaultAssay(tc) <- 'integrated'
tc <- ScaleData(tc, verbose = FALSE)
tc <- RunPCA(tc, verbose = FALSE)
tc <- RunUMAP(tc, dims = 1:50)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
library(dittoSeq)
dittoDimPlot(tc, var = 'orig.ident', size = 0.25)
dittoDimPlot(tc, var = "Phase", size = 0.25)
dittoDimPlot(tc, var = "celltype_fine", size = 0.25)
tc <- Seurat::FindNeighbors(tc, dims = 1:50, reduction = "pca")
tc <- Seurat::FindClusters(tc, resolution = seq(0.5, 1, 0.1))
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 35543
## Number of edges: 1811733
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8738
## Number of communities: 23
## Elapsed time: 12 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 35543
## Number of edges: 1811733
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8641
## Number of communities: 23
## Elapsed time: 12 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 35543
## Number of edges: 1811733
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8552
## Number of communities: 25
## Elapsed time: 12 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 35543
## Number of edges: 1811733
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8473
## Number of communities: 26
## Elapsed time: 12 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 35543
## Number of edges: 1811733
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8401
## Number of communities: 26
## Elapsed time: 12 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 35543
## Number of edges: 1811733
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8324
## Number of communities: 26
## Elapsed time: 13 seconds
pdf("plots/umap_t_cells_cl0.8.pdf", width = 9)
p <- dittoDimPlot(tc, var = 'integrated_snn_res.0.8', size = 0.5, do.label = TRUE)
print(p)
dev.off()
## quartz_off_screen
## 2
print(p)
tc$tcell_annotation <- "other"
DefaultAssay(tc) <- "RNA"
markers <- read.csv("metadata/marker_genes_man_tcell.csv")
mli <- split(markers, markers$tcell_annotation)
for(cellt in names(mli)) {
marker_genes <- mli[[cellt]]$gene
cellt_name <- gsub("[ \\+]", "", cellt)
pdf(paste0("plots/UMAP_tcells_markers_", cellt_name, ".pdf"))
multi_dittoDimPlot(tc, var = marker_genes,
min.color = "grey",
max.color = "purple",
size = 0.25,
order = "increasing")
dev.off()
}
dittoDotPlot(tc, vars = unique(markers$gene), group.by = "integrated_snn_res.0.8")
tc$tcell_annotation[tc$integrated_snn_res.0.8 %in% c(3, 12)] <- "NK cells"
tc$tcell_annotation[tc$integrated_snn_res.0.8 %in% 10] <- "T regulatory cells"
tc$tcell_annotation[tc$integrated_snn_res.0.8 %in% 9] <- "T effector cells"
tc$tcell_annotation[tc$integrated_snn_res.0.8 %in% c(5,8)] <- "Naive T cells"
tc$tcell_annotation[tc$integrated_snn_res.0.8 %in% c(1,7,15,4)] <- "Cytotoxic T cells CD8+"
tc$tcell_annotation[tc$integrated_snn_res.0.8 %in% c(2,13)] <- "Exhausted T cells CD8+"
p <- dittoDimPlot(tc, var = "tcell_annotation", do.label = TRUE, size = 0.25,
labels.size = 2)
print(p)
pdf("plots/umap_manual_annotation_tcells_labeled.pdf", width = 10)
print(p)
dev.off()
## quartz_off_screen
## 2
p <- dittoDimPlot(tc, var = "tcell_annotation", size = 0.25)
print(p)
pdf("plots/umap_manual_annotation_tcells.pdf", width = 10)
print(p)
dev.off()
## quartz_off_screen
## 2
p <- dittoBarPlot(tc, var = "tcell_annotation", group.by = "orig.ident",
x.reorder = match(sample_order$sample,
unique(tc$orig.ident)))
print(p)
pdf("plots/barplot_manual_annotation_tcells.pdf", width = 10)
print(p)
dev.off()
## quartz_off_screen
## 2
library(dittoSeq)
p <- dittoDotPlot(tc, vars = unique(markers$gene),
group.by = "tcell_annotation",
min.color = "lightgrey", max.color = "darkred")
print(p)
pdf("plots/dotplot_tcell_markers.pdf")
print(p)
dev.off()
## quartz_off_screen
## 2
tmark <- read.delim("metadata/tcell_markers.txt", header = F)
cr <- colorRampPalette(c("white", "black", "red"))
p <- dittoHeatmap(tc, unique(c(markers$gene, tmark$V1)),
scaled.to.max = TRUE,
heatmap.colors.max.scaled = cr(25),
annot.by = c("tcell_annotation", "integrated_snn_res.0.8"))
pdf("plots/heatmap_tcell_markers.pdf", height = 10, width = 15)
print(p)
dev.off()
## quartz_off_screen
## 2